Few-particle Green's functions for strongly correlated systems on infinite lattices 



Mona Bcrciu 

Department of Physics & Astronomy, University of British Columbia, 

(Dated: December 9, 2011) 



Vancouver, BC, Canada, V6T IZl 



We show how few-particle Green's functions can be calculated efficiently for models with nearest- 
neighbor hopping, for infinite lattices in any dimension. As an example, for one dimensional spinless 
fermions with both nearest-neighbor and second nearest-neighbor interactions, we investigate the 
ground states for up to 5 fermions. This allows us not only to find the stability region of various 
bound complexes, but also to infer the phase diagram at small but finite concentrations. 
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PACS numbers: 71. 10. Li, 31.15.ac, 71.35.Pq 

Recently, there has been considerable interest in few- 
particle solutions of interacting Hamiltonians. For exam- 
ple, in Ref. [1] it was shown that knowledge of the two- 
and three-body solutions allows for quantitatively ac- 
curate predictions of finite-temperature thermodynamic 
quantities for many-body systems. As another example, 
in the context of atomic and molecular physics, the pre- 
dicted universal three-body Efimov structures [2] have 
now been seen experimentally [■->] , giving new impetus to 
their study and work on various generalizations [4]. 

While the above work is for free space where particles 
have parabolic dispersions, there is equally strong inter- 
est in the lattice version of such few-body problems. For 
example, while stable excitons - bound pairs comprised 
of an electron and a hole - appear in many materials, 
it is less clear when a so-called charged exciton or trion, 
consisting of two holes and one electron or viceversa, is 
stable. That this can happen has been recently demon- 
strated in GaAs quantum wells [5] and in carbon nan- 
otubes [(i]. (Note that trion theory is still mostly based 
on continuous models and variational solutions, e.g. see 
Ref. [7]). Studying bigger bound complexes, for example 
bi-exciton pairs, is the next logical step. 

Few-particle bound states are relevant not only for 
the materials where they appear, but also in the inter- 
pretation of certain spectroscopic data. For instance, 
the role played by bound two-particle states, leading to 
atomic-like multiplet structures in the Auger spectra of 
narrow-band insulating oxides, is well established [^]. At 
low dopings, more complicated complexes may form and 
leave their fingerprints in various spectroscopic features. 
It is therefore useful to be able to study relatively easily 
few-particle solutions on an infinite lattice. 

In this Letter we show that few-particle Green's func- 
tions can be calculated efficiently for strongly correlated 
lattice Hamiltonians in the thermodynamic limit, at least 
so long as the hopping involves only nearest neighbor 
sites. For simplicity and to illustrate the technique and 
its usefulness, we focus here on a one-dimensional (ID) 
model of spinless fermions with nearest-neighbor (nn) 
and next-nearest-neighbor (nnn) interactions. However, 
the method generalizes straightforwardly to higher di- 
mensions, longer (but finite) range interactions, mixtures 



of fermions (including spinful fermions) and/or bosons, 
etc. Such problems are of direct interest either in solid 
state physics, or for cold atoms in optical lattices. 

For two-fermion Green's functions {Nf — 2), our 
method is equivalent to that of Ref. [S] , but is recast in a 
simpler form which allows, in ID, for an analytical solu- 
tion for any finite-range interaction. More importantly, 
it has a simple generalization for Nf > 2. We study cases 
with up to Nf = 5 and show that these suffice not only 
to sort out the stability of few-particle bound states, but 
also to infer the low density phase diagram. 

Consider, then, spinless fermions on a ID chain with 
N ^ oo sites, described by the Hamiltonian: 



n 



h.c. 



H+l 



U2 E nini+2 



where Ci removes a spinless fermion from site i located at 
Ri = ia, and rii = c\ci. Note that this ID Hamiltonian 
is not integrable in the sense of having a Bethe ansatz 
solution. Because our solution is not linked in any way 
to such integrability, it can be generalized to higher di- 
mensions, as already mentioned. To illustrate the main 
idea behind our solution, we discuss in some detail the 
solution for Nf = 2 fermions, after which we generalize 
to Nf > 2. Other possible generalizations, mentioned 
above, arc discussed in the supplementary material [9]. 

Because the Hamiltonian is invariant to translations, 
the total momentum of the pair is a good quantum num- 
ber. As a result, we work with the Nf ~ 2 states: 



\k,n) = — ^ 

Vn 



which describe fermions at a relative distance n > I. 
We define the two-particle Green's functions: 

G{m,n;k,uj) = {k,m\G{uj)\k,n) 

where G(uj) = [ui + iri — TL]^^ with ?/ — > 0+ and wc set 
h ^ 1. From the Lehmann representation: 



G(rn, n; k, w) = 



{k, m\k, a) {k, a\k, n) 
w - E2,a{k) + ir] 



where {\k, a)} are the two-particle eigenstates with total 
momentum fc, 'H|fc,a) = E2^a{k)\k^ a) . Thus, this prop- 
agator allows us to find the Nf = 2 spectrum and also 
to get information about its eigenfunctions. Its Fourier 
transform G{m,n; k,t) (x (fc, m| exp(— n) is the 
amplitude of probability that if initially the two particles 
(with total momentum k) are at a relative distance na, 
they will be at a relative distance ma after time t. 

Matrix elements of the identity 1 = G(aj) {u + irj — H) 
lead to (5„,„i = {lo + ir])G{m, n;k,uj) — (k, m|G(w)'H|fc, n). 
Smcen\k,n) = U{n)\k,n) - f{k) 1) + \k,n+l)], 

where U{n) = U\bn,\ + C^2<^n,2 and /(fc) — 2tcos^, we 
get a simple recurrence relation: 

<5n,m = [uj + iri - U (rt)]G(m, n; fc, w) 

+ f{k) [Gim,n - l;k,uj) + Gim,n + l;k,uj)] (1) 

This is trivial to solve for an infinite chain if one realizes 
that for any m of interest, G{m, n; fc, oj) — )■ as ri — >■ oo. 
This is obvious if to is outside the free two-particle con- 
tinuum where eigenstates, if any, are bound and there- 
fore wavefunctions decay exponentially with n. It is also 
true inside the free two-particle continuum. Even though 
here the wavefunctions are plane- waves, 77 defines an ef- 
fective lifetime r ^ I/77. As such, G(m,n;fc,i) — > if 
na is large compared to the typical distance that parti- 
cles travel within r. Thus, the recurrence relation can be 
solved starting from G{m,Mc + l;k,Ld) ~ for a suffi- 
ciently large cutoff Mc- Of course, the Nf = 2 case can be 
solved analytically exactly (see below). However, the idea 
can be used for Nf > 2 cases, where a numerical solution 
is needed. Noting that the few-particle Green's functions 
become arbitrarily small as a "relative distance" M (to 
be defined below) increases, the recurrence relations can 
be solved propagating the solution from a cutoff Mc to- 
wards small M. Mc is then increased until convergence 
is reached. The effects of rj and Mc on the numerical 
solution arc discussed in the supplementary material. 

First, though, we complete the Nf — 2 discussion, 
which has an analytical solution (for details see [9]; we 
also show there how to deal with a finite-size system 
in this case). At the Brillouin zone (BZ) edge, since 
f{k = 7r/a) = we find: 

G{l,n] -,uj) 



OJ + ir] — Ui' 



as expected since If , 1) is an eigenstate of H with energy 
Ui. For any fca 7^ tt, we find 



G(l,l;fc,w) = 
and for any n > 2 
G(l, n; k, lo) 



Lu + irj — Ui — 



LO + i-q - U2 + z{k,Lo)f{k) 

[z(fc,c^)]"-i/(fc)G(l,l;fc,c.) 
LO + ~ U2 + z{k, Lo)f{k) 




FIG. 1: (color online) A2{k,uj) for U2 = and (a) Ui = 0; (b) 
Ui = —1.5t and (c) Ui = —3t. As the nn attraction is turned 
on, a bound state splits from the free two-particle continuum 
shown in (a). It exists at all momenta if Ui < —2t, but only 
for large momenta if Ui > — 2f. Here 77 — 0.01. 



Values for m > I can be obtained similarly. Here, z{k,Lo) 
is the root of the characteristic equation of this recurrence 
relation, {lo + irj) + f (k) (z -|- i) = 0, for which |z(fc, cj)| < 
1 [!)]. This shows that indeed, G(l, n; , fc, w) — > as n — > 
00. It is also easy to check that inside the free two- 
particle continuum, |a;| < 2/(fc), we have 1 — |2:(fc,a;)| ~ 
?7, so here G decays exponentially only because rj > 0. 

To study the two-particle spectrum, we plot the two- 
particle spectral weight: A2{k,Lo) = — ilmG(l, 1; fc, w) in 
Fig. (1) for U2 = 0, and three values of Ui. By definition, 
A2{k, lo) is finite at energies in the two-particle spectrum, 
and its value is related to the probability to find the 
fcrmions as nn in that eigenstate. If Ui = 0, A2{k,Lo) 
is finite in the free two-particle continuum, ranging from 
—At to 4t if fc = 0, while at fc = n/a only w = is an 
eigenstate, hence the (5-function (Lorentzian) seen here. 
As an attractive Ui is turned on, the fc = 7r/a peak tracks 
Ui, and a bound state is pulled below the continuum at 
nearby k values. For Ui > —2t, this bound state exists 
only near the BZ edge, while near the F point the weak 
attraction shifts spectral weight to the bottom of the two- 
particle continuum but is not enough to push a discrete 
state below it. For Ui < ~2t, the bound state becomes 
the low-energy state at all fc. This shows that, for certain 
ranges of parameters, bound pairs are only stable in some 
regions of the BZ, which moreover are not necessarily 
near fc = 0. It would be interesting to investigate their 
effects on various response functions. 

However, hereafter we focus on the fc = ground- 
-1 state (GS). Fig. 2a shows whether in the GS the pair 
is bound or not, for Ui < and U2 > 0. (Note that 
such interactions, attractive at short-range and repulsive 
at longer-range, appear in systems with highly polariz- 
able ions [10]). For Ui < —At a bound pair is always 
stable; even if it had infinite mass, a nn pair of energy 
Ui is below the minimum energy of two free fermions. 
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FIG. 2; (color online) Stability diagram for (a) Nf = 2, and 
(b) Nf = 3 fermion systems, indicating the nature of the GS. 
The inset in (a) shows that bound pairs are always stable if 
Ui < —At. The dashed line is a perturbational prediction. 
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FIG. 3: (color online) Nj — 5 GS energy (green circles, label 
"5") vs. U2 for (a) Ui = -3.5t, and (b) Ui = -2.5t. Other 
lines show the lowest energies of various complexes, and ar- 
rows indicate dissociations (see text for details). 



of —At. Of course, the kinetic energy of the pair further 
enhances its stabihty region. The full line in the inset 
shows a perturbational estimate for i <C |C/i — L/2I [■■)]. 

This A^^ = 2 stability diagram, however, has no pre- 
dictive power for what happens if more fermions are in 
the system. For example, if Nf — 3, wc expect regions 
where the GS consists of 3 fermions, of a bound pair 
plus a fermion, or of a bound "trion" . To identify these 
regions we study Nf = 3 Green's functions, by direct 
generalization of the Nf = 2 approach. Briefly, for any 
ni > 1,712 > 1, we define three-particle states: 

^ i 

and three-particle Green's functions: 

G{mi,m2;ni,n2;k,uj) = (fc, mi, m2|G(w)|fc, ni, 712). 

Recurrence relations for these propagators arc generated 
just as for the Nf — 2 case. If we define a "relative dis- 
tance" M = rii -f 77,2, hopping of the outside fermions 
will link Green's functions with a given M to those with 
A/±l. If the central fermion hops, one of the 77i, 712 values 
increases by one and the other decreases by one, there- 
fore M remains the same. Thus, the equation of motion 
links Green's functions with consecutive M — 1 , M, M + 1 
values, leading to recurrence relations that can be solved 
in terms of continued fractions of matrices, if we use the 
insight that propagators vanish as M ^ 00. Generaliza- 
tion to larger Nf values is now straightforward [')]. 

In higher dimension, we need to combine the "relative 
distance" with the "Manhattan distance" [11]. For exam- 
ple, in 2D for Nf = 3, wc associate the plane-wave with 



the coordinates ix and iy of the "central" particle for 
that axis. The other particles' coordinates are ix — ni^x, 
ix+n2,x, respectively iy-ni^y, iy + n2,y, where m^a > 0, 
i = 1,2, a = x,y. If wc choose M = X^i a then 
nn hopping links together only Green's functions with 
M — 1, M, M+1. Thus, m particles in 2D is computa- 
tionally similar to 27n — 1 particles in ID. In both cases, 
2{m — 1) positive integers specify the relative positions, 
and M is their sum. The key observation is that the 
equations of motion still group into recurrence relations 
linking only quantities with M — 1,M,M + 1, allowing 
for an efficient solution (for more details, see [9]). 

To study the spectrum of the Nf = 3, ID system, we 
plot ^3(^,0;) = — ilmG(l, 1; 1, 1; fc, w). This must have 
finite spectral weight for lu > i?2,GS— 2t, corresponding to 
a continuum of states describing a fermion far away from 
a pair. (If i?2,GS = ~4i, this continuum starts at —6t and 
describes 3 free fermions) . If the continuum is the lowest 
spectral feature, then the GS is either a pair+fcrmion or 
three fermions, mirroring the Nf — 2 situation. However, 
if a discrete state appears below this continuum, then the 
GS is a stable bound trion [9]. The stability diagram is 
plotted in Fig. 2b and shows a region where trions are 
stable, at large attractive \J\ and weak repulsive U2- This 
is expected since binding a 3rd fermion to a stable pair 
lowers its energy by roughly Ui-\-U2, while a free fermion 
can lower the total energy by at most —2t. 

The fact that stable trions are found for Nf = 3 does 
not, however, guarantee that they appear at finite con- 
centrations. Just as the pair-l-fermion is unstable to trion 
formation, trions may be unstable to bigger bound com- 
plexes, if more particles are present. Indeed, a study of 
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FIG. 4: (color online) Phase diagram at small concentrations. 
The GS consists either of unbound fermions, or of pairs, or it 
phase separates into fermion rich and fermion poor regions. 
The dashed line is an estimate for phase separation (see text). 



cases with Nf ~ 4 and 5 fermions proves that trions are 
actually unstable. This is shown in Fig. 5(a) where we 
plot the energy of the iV/ — 5 GS vs. U2 (hne marked 
"5") at U\ = — 3.5t. The other lines show energies where 
a continuum could appear, eg. E2+2+1 = 2£'2, GS — 2< is 
the lowest energy of two pairs plus a fermion, £'2+3 ~ 
E2,GS + E^^as is the lowest energy for a pair plus a 
trion, etc. The arrows indicate various dissociations. Ar- 
row 1 shows when a pair becomes more stable than 2 
fermions (i?2+i+i+i < ^'i+i+i+i+i), while arrow 2 shows 
when a trion becomes more stable than a pair+fermion 
(£■3+1+1 < £2+1+1+1), see Figs. 2a, b. A trion+fermion 
is unstable to either two pairs (at larger U2) or a 4- 
fermion bound complex (smaller 1/2)- The boundary be- 
tween the two is marked by arrow 3 (£4+1 = £2+2+1)- 
But 4-fermion bound states are not stable either, since 
£4+1 < TOzn(£2+3,£5) (arrow 4 marks where the 5- 
fermion bound complex breaks into a pair+trion). Be- 
low it, £5 is indeed in good agreement with the pertur- 
bational estimate for the energy of a 5-bound complex 

£5,5 = 4c/i + 3U2 + 2ty{Ui + e/Ui - 2t^/{Ui + 1/2)), 

shown by the dashed hne indexed "5, bound". 

What happens as Nf increases becomes clear if we re- 
alize that arrows 3 and 4 point to essentially the same 
U2 value. If more fermions are added, below this U2 
we expect a bigger and bigger bound complex - in other 
words, phase separation occurs and the system splits into 
a fermion rich and a fermion poor region. Above this, a 
gas of pairs is stable (plus one trion, if Nf is odd) . That 
this inference is correct is verified by the following argu- 
ment. This critical value should be given by the condition 
that adding two more particles to a fermion rich region 
(which changes energy by about 2Ui + 2U2, because of 
extra interactions) should be energetically favorable to 
having a bound pair far away. From 2U2 + 2Ui < £2,05 
we find U2 = 1.29t if Ui = —3.5t, in good agreement with 
the value U2 = I St pointed to by arrows 3 and 4. 

Thus, based on these few-particles results, we can infer 



the phase diagram of this model at small concentrations, 
shown in Fig. 6. The dashed line shows the estimate 
discussed above, accurate for large Ui,U2 (at smaller Ui, 
t comes into play since the extra fermions need not be 
fully localized at the edge of the fermion rich region). 
If Ui > —2.6t, the transition is from phase separation 
to unbound fermions as U2 increases. This is shown, for 
Ui = — 2.5t, in Fig. 5b: here each bigger complex is more 
stable than any smaller ones, if U2 < 0.63t (arrow). 

While we are not aware of numerical studies of this 
model, the good agreement with various asymptotic esti- 
mates as well as with known results for spin-^ Hamilto- 
nians [9], supports the accuracy of our results. This work 
shows that even such a simple model has a rich behavior 
that can be uncovered with this method. 

To summarize, we have shown how to calculate few- 
particle Green's functions on an infinite ID chain. The 
information obtained from them sheds light on the sta- 
bility of few-particle bound states. It also illustrates the 
dangers of an insufficient analysis ~ if we stopped at 
Nf = 3, we would conclude that trions are stable in a 
large region of the parameter space, in this model. Anal- 
ysis for larger Nf shows that addition of more particles 
leads to instability of trions, and furthermore allows us 
to find the phase diagram for small concentrations. 

Although these results are for a ID model, as discussed 
above this method generalizes to higher-D if the hopping 
is nearest-neighbor only. This opens the way to study 
the stability of trions and bi-excitons in realistic lattice 
models. Such work is currently under way. 
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SUPPORTING MATERIAL 
Further support for this phase diagram 

As discussed in the main text, the phase diagram we 
derived from the few-particle Green's functions is sup- 
ported by the asymptotic hues shown in Figs. 2a (inset) 
and Fig. 4. We are not aware of any numerical studies of 
this model that could be used for direct comparison, al- 
though it bears mentioning that once is large enough 
that convergence has been achieved (see discussion and 
examples below), these results are exact - there is no 
approximation involved in obtaining them. 

However, additional support for these results can be 
obtained from studies of spin— i Hamiltonians. Through 
a Jordan- Wigner transformation [1], the Hamiltonian for 
spinless fermions studied here can be mapped into: 

n = Y, [-2t {sfsf^, + si^i) + c/iSf -I- u^stst^^ 



them to be the shortest possible "arc" , n = min(n, N—n). 
Using e**^^" = 1, it follows that for any 1 < n < Nmax 



\k,N- 



Ifc, 



The equations of motion for the two-particle Green's 



functions were derived in the main text. For m : 
using the notation G(l,n; k,uj) — )■ a„, they are: 

(w + 17] - Ui)ai + f{k)a2 = 1 

{uj + iri- U2)a2 + f{k){ai + 03) = 

{uj + iri)an + /(fc)(a„_i + a„+i) = 

for 3 < n < Nmaxi and finally 
(w + ir])aN^^^^ 



fik) 



1 - e- 



■ika 



aw 



1, and 

(2) 
(3) 
(4) 

(5) 



A numerical solution is now trivial. 

However, we can do better. First, note that for ka = ir 
we have f(k) = 0, and so the solution is: 



The case Ui < 0,U2 > corresponds to ferromagnetic 
nearest-neighbour Ising interaction and frustrating anti- 
ferromagnetic next nearest-neighbour Ising interaction. 

For 1/2 = 0, the phase diagram of this model is well- 
known: it consists of a ferromagnet for \Ui\ > 2t and 
a Luttinger liquid otherwise [2]. The ferromagnet corre- 
sponds to phase-separation in the fermion language, since 
the magnetization is linked to the density of fermions. 
The unpaired fermions phase is reasonably linked to the 
Luttinger liquid, so this line of the phase diagram agrees 
with other known results. 

While we could not find a study of the spin model with 
[72 ^ 0, it is expected that an increase in U2 frustrates the 
ferromagnetic phase and eventually makes it unstable. 
We believe that the pairs phase is a charge 2 Luttinger 
liquid, a sort of ID version of a p-wave superconductor. 
Such phases are fairly well known [2, 3]. 



Details for the Nf = 2 solution 

As discussed, because the Hamiltonian is invariant to 
translations it is convenient to use the two-fermion states: 



G(l,?i; -,uj) 



The method can be trivially extended to systems with 
disorder because using these translational states is not 
an essential ingredient of the method. 

For a finite-size chain with N sites, in order to not 
double count the states, we must restrict 1 < n < Nmax, 
where N,nax = if iV is even, respectively Nmax = ^^y^ 
if N is odd. In other words, when considering the two 
fermions on the closed ring, we take the distance between 



+ 4?/ - Ui 



This is expected, since 1) is an eigenstate of the full 
Hamiltonian with energy Ui , as can be easily checked. 

For any ka 7^ tt, the equations for 3 < 71 < Nmax — 1 
can be solved analytically. Their general solution is: 



where z and are the roots of the characteristic equa- 
tion: (w + irj) + f{k) {z + ^) = 0. We choose z = z{k, uj) 
to be the root for which \z\ < 1. This is always uniquely 
defined for any > 0, since the product of the two roots 



fik) 



fik) 



is 1. The solution is now immediate. The last equation 
for n = Nmax fixes the ratio of (3 /a, the solution is then 
propagated down to n = 2, and one uses the first two 
equations to find ai and a, completing the solution. 

These results have a very straightforward physical in- 
terpretation. Note that the general recurrence equation 
for 3 < n < Nmax is the same that describes two free 
fermions. For a given k, the two-fermion continuum 
spans the energies {ek-q + eq}q = [— 4i cos 4t cos 
where = — 2tcos(A:a) is the free particle energy. 

This explains why for energies outside this range |w| > 
2f{k) = 4tcos^, we find the two roots z,l/z to be 
real (when 77 — >■ 0), and such that \z\ -> as |a;| — ?> 00. 
This shows that here the two-particle Green's function 
decays exponentially with n, as expected since there are 
no free two-particles eigenstates at these energies. The 



exponentially increasing 



part is a finite size 
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effect: on a finite ciiain, the inter-particle distance even- 
tually decreases as n N. In the limit iV — >■ oo, this 
contribution must vanish, in other words in the thermo- 
dynamic limit we must have a„ = zUn^i = az"" as a 
purely exponentially decreasing function with distance. 
This is consistent with the fact that the solution must be 
insensitive to how we end the recurrence relation (what 
is the boundary condition) as TV — >■ oo. 

For certain values of Ui,U2, new eigenstates may ap- 
pear outside the free two-particle continuum. These de- 
scribe bound pairs, therefore we expect their wavefunc- 
tions (and the associated Green's functions) to decay ex- 
ponentially with the distance n between the particles. 
This is fully consistent with the previous discussion. In 
particular, since |z| — )■ as w — > — oo, it follows that 
lower-energy pairs are more tightly bound. 

On the other hand, inside the free two-particle con- 
tinuum, |a;| < 2/(fc), wc find |z| — > 1 as 77 — > 0. As 
a result, the two contributions to a„ become oscillatory 
functions, underlying the fact that there are freely prop- 
agating two-particle eigenstates at these energies. The 
interaction will be responsible for scattering leading to 
phase-shifts, however at these energies the two particles 
can propagate arbitrarily far from each other, if ?/ = 0. 

Using a finite rj (which wc arc forced to do, for numer- 
ical reasons), results in 1 — |z| ^ 77, in other words there 
is slow exponential decay at these energies as well, but 
now controlled by 77. Physically, this is because the finite 
77 introduces a "life-time" for these particles. As a result, 
if the distances n,N — n are very large compared to the 
typical relative distance the two free particles can explore 
in their lifetime t ^ the probability for the pair to be 



at such distances decreases exponentially, and so do the 
Green's functions. It follows that in the thermodynamic 
limit, we can again take a„ = az". Then, the solution in 
the thermodynamic limit is trivial since we only need to 
solve the equations for ai and 02 using 03 = z(fc,a;)a2. 
The solution is listed in the main text. Generalization 
to longer (but finite) range interactions is trivial, as is 
finding the solution for other m values. 



Nf = 3 in the thermodynamic limit 

Based on the arguments discussed above, in the limit 
of an infinite chain we expect the Green's functions to de- 
cay exponentially at all energies, with an exponent con- 
trolled by 77 inside the three-particle continuum, and by 
the inverse of the distance between w and the continuum's 
band-edge, for energies outside the continuum. 

We use three-particle states of total momentum k: 

i 

where rii > 1, 712 > 1. Since we take iV — 00, the states 
with ni ^ Nmax,n2 Nmax bccomc irrelevant and we 
do need to worry about properly counting them; because 
of the finite lifetime, even inside the continuum particles 
cannot travel that far from each other. For a finite chain, 
however, proper counting is important (see below). 

The equations of motion for the three-particle Green's 
functions defined in the main text, are: 



Sni,miSn2,7n2 = {uJ + IT] - f7„i ,„2 )G(7ni , 7712 ; rii , 712 ; k,uj) 

+t [G(777i, 7712; rii — 1, 772; fc, w) + G{mi, 70,2; n-i + 1, 77,2; fc, w) 
-|-e*''''^G'(7?7i, m2; ni — 1, 712 + 1; k, uj) + e~''^°G(TOi, TO2; ni + 1, 712 — 1; k,uj) 

+ G(777l, 7712 ; 771, 712 + l;k,Uj) + G(7?7l , 7772 ; 77l , 772 - '^;k,Uj)]. 

I 



Here, J7„i,„2 = C^l(<5ni,l + l^na.l) + t^2('^ni,2 + <5„2,2 + 

Sni+n2,2) is the interaction energy when the three parti- 
cles are at relative distances ni, 77,2 from each other. The 
remaining terms describe the effect of hopping on the 
771,772) state. The first two terms come from the hop- 
ping of the left- most particle, which changes ni . The next 
two terms come from the hopping of the central particle, 
which keeps ni -f 77.2 constant, and the last two terms are 
from the hopping of the rightmost particle, which changes 

772. If 77l = 1 then G(777l, 7772; TT-l — 1, 772, k,Uj) =0 sluCC 

this hopping process is not allowed for spinless fermions, 
and similarly for 77,2 . For N ^ 00 we need not worry what 



happens as 7ii,7i2 ~ Nmax, since the Green's functions 
vanish before the particles go so far from each other. 

Suppose we are interested in 7771 = 777,2 = 1 and use the 
shorthand notation 0(711,772) = G(l, 1; tii, 712; fc, w). The 
resulting infinite (in the thermodynamic limit) system of 
coupled recurrence relations can be solved as follows. We 
define the vectors 

/a(l,A/-l)\ 
Vm^ a(2,M-2) 

Va(M- 1,1)/ 
which collect all the Green's functions with the same "rel- 
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ative distance" M = ni + n2- Its dimension is M — 1, 
although for special values of k there are further symme- 
tries that can lower it. For example, at fc = we have 
a(n, M — n) = a{M — n, n) and the dimension is halved. 

The special property of Hamiltonians with only 
nearest-neighbor hopping is that the resulting equations 
of motion only link three vectors with consecutive rela- 
tive distances. In other words, for any M > 3, we can 
recast the recurrence equations as: 

7a/ Va/ = aAiVM-i + I3aiVai+i (6) 

where aM , Pm , 1m are very sparse matrices whose matrix 
elements are simple functions of fc,a; that can easily be 
read off the equations of motion. Because we know that 
all Green's functions must vanish in the limit M — >■ oo, 
the solution of this recurrence equation is given by: 

Va/ — AmVm-i 

where the matrices Am are given by continued fractions: 

Am = [iM - PmAm+iI^^olm 

and can be calculated starting with Am^+i = at a 
sufficiently large cutoff Mc- 

Once all these matrices are known, and in particular 
A3 which links a(l,2) and a(2, 1) to a(l,l), we can use 
the equation of motion with rii = 71,2 = 1 to find: 



' ' ' CJ + i77-2C/i -[/2 + t[A3|l,l+^3|2,l]' 

from which we can then get all the other propagators. 

To illustrate the effect of the numerical parameters 77 
and Mc, we analyze the three-particle spectral weight: 

A3H = --ImG(l,l;l,l;fc = 0,c^) (7) 

TT 

This is finite for all energies w in the fc = 0, Nf = 3 
spectrum, and its weight gives the probability of having 
the three fermions located on three consecutive sites. 

Just like for Nf = 2, whether a bound trion is the 
ground-state or not is determined by whether a discrete 
Lorentzian appears below the continuum, or not. The 
continuum starts at i?2,GS ~ 2t. where -E2,gs is the GS 
energy of the Nf ~ 2 case. If the parameters arc such 
that bound pairs are not stable, then £2^3 = and 
the three-particle continuum starts at — 6i. However, if a 
bound pair is stable, then the continuum moves to lower 
energies, and consists of states where a free particle scat- 
ters off a bound pair (other higher-energy features are 
also present, but not of interest for our analysis). 

In Fig. (5), we show ^3(0;) for Ui = —3t, U2 = and 
three sets of parameters 77, Nc- The dashed vertical line 
shows the expected on-set of the continuum, at E2,GS — 2t 
[for these parameters, E2,gs ~ — 4.33t, see Fig. Ic in the 
main text]. Clearly, A^iuj) shows a continuum starting 




FIG. 5: ^43(01) vs uj for Ui = —Zt, U2 = 0, for various values 
of the broadening 77 and the cutoff Mc- The dashed line shows 
the expected continuum onset at E2.GS — 2t. 



at this energy, but there is also a Lorentzian peak below 
it, indicating a stable trion for these parameters. 

Note that the spectral weight in the continuum de- 
pends on the specific broadening rj and cutoff Mc used. 
This dependence can be understood easily. Afc is the cut- 
off at which we set the Green's functions to zero. Physi- 
cally, this is equivalent with adding an effective "interac- 
tion" which becomes infinite if the total relative distance 
between particles n + m > Mc, and is zero otherwise. 
As is the case for any system in a "box" , we expect the 
continuum to be replaced by a set of discrete levels, with 
a spacing 5E ~ 1/A/c- These states, however, are broad- 
ened by 77. This explains why the first curve is much 
smoother than the second one, even though they have 
the same Mc- On the other hand, the third curve has 
Mc increased by a factor of two, and indeed there are 
roughy twice as many oscillations marking the discrete 
peaks. For any value of Mc, the curve becomes smooth if 
77 is large enough so that 5E ^ -q. As already discussed, 
physically this means that the lifetime t 1/r; is so short 
that the particles cannot travel up to the boundaries of 
this potential "box" defined by Mc- 

Below the continuum, the spectral weight is insensitive 
to Mc, because the states that appear here (if any) are 
bound well inside this "box" . The broadening 77 is still 
reflected in the shape of the Lorentzian: although not 
shown entirely in Fig. 5, the peak for the smaller 77 is 5 
times narrower and taller, as expected. This insensitivity 
to Mc is very convenient, because it means that one can 
get very good estimates for the energy of strongly bound 
states using rather small Mc values. Of course, if the 
binding energy is very small, then one has to increase 
Mc until convergence is achieved. 
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FIG. 6: CPU time for one frequency, vs. cutoff Mc, for sys- 
tems with Nf = 3, 4, 5. 

Higher Nf in the thermodynamic Umit 

It should now be apparent that the method generahzes 
straightforwardly to any case with an odd number Nf of 
particles. We choose the reference location i as being 
that of the central particle, and index states in terms of 
the absolute values of the relative distances of all other 
particles ni, n2, njvj-i with respect to the central par- 
ticle. In the equation of motion, hopping of any of these 
other particles will increase or decrease its own rii by 1, 
so the "relative distance" M = 'Yl!i=i ^ varies by 1. If 
the central particle hops to the right, for example, this 
increases by 1 all the distances to all particles to its left, 
and decrease by 1 all distances to all particles to its right. 
Since there are equal numbers of particles to the left and 
to the right of the central particle, then M is unchanged. 

For an even Nf, we choose as the reference particle 
either of the two central particles. In this case, the "rel- 
ative distance" is changed by 1 when any of the particles 
hop, including the "central" one. 

In either case, the recurrence equation can still be cast 
in the general form of Eq. (6) and can be solved by 
similar means. Of course, the larger Nf is, the larger 
is the dimension of Vm, so eventually one runs out of 
computational power to calculate the continued fractions 
numerically. This is the factor that limits what values of 
Nf can be considered. 

Fig. 6 gives the real time to calculate the spectral 
weight at one frequency on a 4-core CPU, for systems 
with Nf = 3, 4 and 5 fermions. As expected, the run 
times increase quite fast with both Nf and the chosen 
cutoff Mc- Note that for Nf = 3,4 we showed data for 
Mj, much larger than what is needed to achieve conver- 
gence, simply because for smaller values the CPU time 
becomes independent of M^., showing that it is deter- 
mined by other tasks, not by the computation of the con- 



tinued fractions which is the most time-consuming part 
at large M^- In fact, Mc = 50 sufficed to achieve con- 
vergence even in the most delicate cases discussed in the 
main article, namely where a bound state is very close 
to a continuum [i.e. near a dissociation process). As 
already mentioned, here one needs to use a small r] and 
therefore a larger Mc to be able to separate such close- 
by features. If the energy of the bound complex is well 
below the continuum, on the other hand, much smaller 
Mc suffices and calculations are much faster. 

In practice, it is useful to first use a fairly small Mc 
to quickly scan a large range of energies to see where 
the main features are. Of course, one particularly useful 
characteristic of this calculation is that the spectrum for 
a given number Nf of fermions must have one or more 
continua at energies determined by the spectra with fewer 
fermions, which arc known. This gives not only a chance 
to validate the computation, but also a very useful in- 
dication of where features are expected in the spectral 
weight. If a A^/-bound complex is stable, its correspond- 
ing Lorentz peak is below the lowest-energy continuum 
and must be found by searching for a peak in the spectral 
weight in the infinite range of energies lying below this 
lowest continuum. Even in this case, one may use pertur- 
bation theory as a first guide to where the peak may be. 
As an example, see line "5, bound" in Fig. 3a which pro- 
vides an estimate for the energy of the 5-fermion bound 
complex. To zero order, its energy is E^^b = 4C/i + 3J72, 
since a configuration with 5 fermions occupying consec- 
utive sites has 4 nn and 3 nnn pairs. If Ui,U2 are com- 
parable to t, then one can use perturbation to allow the 
end fermions to hop one site on and off the end of the 
complex; this further lowers the energy by 2P/Ui (note 
that Ui < 0); since the resulting configuration only has 
3 nn and 3 nnn pairs. If needed, 2nd and higher order 
corrections can be obtained by including further possible 
configurations, in a standard fashion. Using such guid- 
ance plus low-A/c scan of a large range of energies, the 
rough position of the peak can be found efficiently, af- 
ter which Mc is increased until the energy of the peak is 
converged to the desired precision. Note that since this 
peak is a Lorentzian with a known broadening 77, as few 
as 2 points close to its maximum suffice to extract its 
maximum and its weight from fitting. 

This is why even though for Nf = 5 and Mc = 50 it 
takes ~ 20s to calculate the spectral weight at one fre- 
quency, one can actually identify the GS energy with high 
precision within very few minutes, for a given value of the 
parameters. This is also why we are confident that this 
type of calculation can be successfully extended to larger 
Nf, especially if clusters with more than 4 CPUs are 
used to further speed up the computation. We stopped 
at Nf = 5 here simply because this value was sufficient 
to deduce the phase diagram for this model. 
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Mix of two different kinds of fermions 

We now briefly discuss the generalization to a mix of two different types of spinless fermions, still on a ID infinite 
chain. Let Ci and di be their corresponding annihilation operators. For spin-i fermions, one can take Ui^^ = c^; Ui^i = 
di. The model Hamiltonian we consider is a direct generalization that used in the main text: 

H = -tc ^(cjci+j + h.c.) + Ui,c'^ nc^^Tici+i + U2.C ^ ric,inc,i+2 

i i i 

- td '^{dld^+i + h.c.) + Ui^d ^ ndsnd,i+i + U2A ^ nd,ind,i+2 

i i i 

+ Uo^ nc.iUd.i + Ui,m ^ {nd.inci+i + h.c.) + U2,m ^ {nd.inc,i+2 + h.c.) . (8) 

i i i 

I 



In other words, each species has both nn and nun inter- 
actions, while the mixed interactions are on-site, nn and 
nnn. Of course, if the two spinless fermions correspond 
to different spin-projections of the same spinful fermion, 
and the interactions are spin-independent, then one ex- 
pects Ui,c = Ui,d = Ui^m etc. 

Consider first a mixed pair with a total momentum k. 
One expects to be able to factorize the two-particle states 
into analogs of "singlet" and "triplet" (with m = 0) 
states, which should not mix with one another through 
hopping, so that each should have its own set of recur- 
rence relations. This is indeed true, however at finite 
momentum these symmetries get mixed and identifying 
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the proper states requires a bit of work. The solution is 
as follows. We define: 

t{k) = tj'^ + tdc-'^ = T{k)e'^'' 

where T{k) ~ \/^? + *1 + Ztc^d cos(fca) and cos(/)fc = 
(tc + id)cos^/r(fc). Then, let: 

|A:,.s,0) = ^ V 
and for any n > 1, 



rfjcl+j |0> 



(9) 
(10) 



The "s" and "t" labels are because at A: = 0, and if these 
are spin-up and spin-down fermions, these states describe 
the usual singlet and triplet combinations. 



It is now easy to check that the recurrence relations 
do not mix "s" and "t" states together, and each set can 
be solved similarly to that for the spinless N f = 2 case. 
This is a double bonus. First, because it keeps the re- 
currence equations simpler, which makes the calculation 
more efficient. More importantly, one can figure out the 
symmetry of the bound states that form, based on which 
Green's functions exhibit poles at those energies. 



More mixed fermions 

Generalization to more particles follows closely. The 
main ingredient, namely that a relative distance can be 
defined, and that hopping only varies it by at most 1, 
stays the same. The complication is that now parti- 
cles of unlike type can pass by each other, so for ex- 
ample for a three-particle calculation with relative dis- 
tance M = n -I- m, one generally has to include all states 



like X;.e 



ifc-Ri „t 



Ad] 



i — ni i+m 
' ^' '-i+m 



|0>, 



|0) and 



|0). Thus, the dimension of the vec- 
tors with a given relative distance is bigger than if the 
particles were identical. Again, careful consideration of 
symmetries (especially at fc = 0) lowers the dimension 
and make the calculation more efficient, besides provid- 
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ing information on the symmetry of the bound states. 
Bosons 

The calculation can be carried over to bosons trivially. 
The main difference is that one can place any number 
of bosons on the same site, so there are additional re- 
currence equations describing states with shorter relative 
distances than possible for fermions. The nn hopping in- 
sures the same general structure of the recurrence equa- 
tions, and in fact for the terms where there are no mul- 
tiple bosons at the same site, the equations are identical 
with those for fermions in similar configurations. 

Higher dimensions 

To obtain this, one has to combine together the idea of 
a "relative distance", described above for multiple par- 
ticles, with that of a "Manhattan distance" which we 
introduced in Ref. 4 to show how to calculate single- 
particle Green's functions in higher dimensions. For ex- 
ample, for two particles on a 2D square lattice, one needs 
two integers n = (nx,ny) to characterize the relative 
distance between the two particles in states of the form 
|k, n) ^ i e*'''^'cjci_,_n|0). Some restrictions apply 
to the allowed values of Ux , Uy so that double count- 
ing is avoided. For example, we can choose rix > 0; if 
Ux ~ then only rij, > is needed, while if rix > 0, ny 
can take both positive and negative values [5]. Nearest 
neighbor hopping will link the Green's function for this 
ket to the ones corresponding to kets with (rix ± l,ny) 
and {nx,ny ± 1). As a result, here we should choose 
M = \nx\ + \ny\ as a. sum of the relative distances pro- 
jected along all the axes - this is a "Manhattan distance" 
characterizing the relative distance between the two par- 
ticles in this state. Nearest-neighbor hopping then pre- 
serves the general structure of linking Green's functions 
with a given M to only others with M±l, and the general 
approach of rewriting the equations of motions in terms 
of continued fractions carries over. In fact, this prob- 
lem is very similar in structure to that of 3 fermions in 
ID, where we also need two integers ni,n2 to character- 
ize each internal arrangement, and where M = ni + n2. 
The main difference is that in 2D, Uy could be negative 
as well, in other words there are roughly twice as many 
states with a given M then for the 3 particles in ID. This 
suggests a corresponding increase in the computational 
time. In reality, even this increase can be eliminated, 
if one explicitly works with pairs of s-wave or d-wave 
symmetry. In the former case, the Green's functions cor- 
responding to a given Ux and zLuy are equal, while in the 
latter case, they have equal magnitude but opposite sign. 
In either case, the actual number of unknowns is halved. 
As a result, the 2D calculation for cither an s-wave or d- 



wave pair is basically equivalent to a ID calculation for 
3 particles. 

The generalization to 3 particles is described briefly 
in the main article. In this case, we need 4 integers to 
describe the relative positioning with respect to the "cen- 
tral" particle (note that different particles can play this 
role, for different axes). In ID, 4 integers are needed to 
describe the relative arrangements of 5 particles. In gen- 
eral, m particles in 2D require 2(m — 1) integers to spec- 
ify relative positioning from a "central" particle, and as 
such, this calculation maps onto a calculation with 2m— 1 
particles in ID. Just as discussed above for to = 2, at first 
sight there are more states with the same total M in the 
2D problem then in the ID one, because some of the 2D 
integers could be negative while the ID ones are all pos- 
itive. However, if symmetries are properly enforced this 
overall multiplication factor can be removed, at least at 
k ^ 0. As a result, one finds not only the spectrum but 
also the symmetry of the bound complex (if, indeed, such 
a bound complex is stable). 

This is why in terms of running times for higher-D, a 
reasonable estimate can be obtained based on running 
times in ID. For example, the 2D, Nf = 2 case discussed 
above is roughly equivalent to a ID, Nf = 3 computation, 
since in both cases the configurations are characterized 
by 2 integers. It follows that stability for trions {Nf = 3) 
and bi-excitons {Nf = 4) can be studied very easily in 
2D, since it involves problems similar to Nf = 5, respec- 
tively Nf = 7 in ID. In 3D, study of trions would be 
equivalent to investigating Nf = 7 in ID, which can cer- 
tainly be accomplished in a reasonable time on a regular 
desktop with a 4-core CPU. To study bi-excitons in 3D 
(equivalent to Nf ~ 10 in ID), it may be needed to use a 
bigger cluster to lower the computation time. Of course, 
if the bound complex is stable, this will be confirmed by 
a small Mc run, in an efficient fashion. 

Longer- range hopping and/or finite chains 

In these cases, it is impossible to recast the equations 
of motion in terms of recurrence equations for consecu- 
tive vectors Vm-i, Vm, Vm+i- For longer range hopping, 
this is because the relative distance will be changed by 
at least ±2 for second-nearest neighbor hopping. Even 
for nearest-neighbor only hopping, in a finite system this 
simple rule is broken for states with particles separated 
by maximum allowed distances. The only exception is for 

= 2 on a chain, where as discussed, the hopping from 
n = Nmax to Nmax -l- 1 is actually mapped into hopping 
to Njnax — 1, up to a phase factor. In all other cases, the 
hopping out of these states with maximally allowed rela- 
tive distances will map into states which can have quite 
different M values. 

One may still obtain a solution for such cases, by solv- 
ing all of them together as a linear system, instead of 
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factorizing them into a recursive set based on their M 
values. While the dimension of this linear system is much 
bigger, the matrix is extremely sparse and can be dealt 
with efficiently by various known algorithms. 

For an infinite system with longer range hopping, if 
one is interested in bound states, then one can set a quite 
small cutoff resulting in a reasonable computational task. 
For a finite-size chain, one needs some physical intuition 
to decide what states (Green's functions) can be removed 
from the calculation, i.e. how to define a "cutoff". We 
are currently investigating such problems. 
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